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We investigate the use of asymptotically null slices combined with stretching or compactification 
of the radial coordinate for the numerical simulation of asymptotically flat spacetimes. We consider 
a 1-parameter family of coordinates characterised by the asymptotic relation r ~ between 

the physical radius R and coordinate radius r, and the asymptotic relation K ^ for the 

extrinsic curvature of the slices. These slices are asymptotically null in the sense that their Lorentz 
factor relative to stationary observers diverges as T ~ R^^^. While 1 < n < 2 slices intersect 
0 < n < 1 slices end at We carry out numerical tests with the spherical wave equation on 
Minkowski and Schwarzschild spacetime. Simulations using our coordinates with 0 < n < 2 achieve 
higher accuracy at lower computational cost in following outgoing waves to very large radius than 
using standard n = 0 slices without compactification. Power-law tails in Schwarzschild are also 
correctly represented. 


I. INTRODUCTION 

Numerical simulations of astrophysical events in gen¬ 
eral relativity are carried out on a numerical domain that 
comprises a central strong field zone surrounded by a 
much larger outer zone which contains no matter, where 
the fields are weak, and in which gravitational waves are 
mostly outgoing. 

One reason for making the numerical domain large is 
mathematical. Boundary conditions (BGs) for the con¬ 
tinuum initial-boundary value problem which are both 
compatible with the Einstein constraints and make the 
problem well-posed have been suggested only recently 0 
and are still untested numerically. BGs which are con¬ 
sistent with the constraints (but not proven to be well- 
posed) Bhave been used successfully in 3D testbed sim¬ 
ulations!^ 0 . Large simulations of astrophysically in¬ 
teresting scenarios have so far been using BGs at finite 
radius (with the exception of j^, see below) that are in¬ 
compatible with the constraints, so that finite constraint 
violations propagate in from the boundary even in the 
limit of infinite resolution. 

A second reason is physical. The central object should 
be isolated with no incoming gravitational radiation. 
This cannot be achieved with purely local BGs. For the 
same reason gravitational radiation cannot be read off 
reliably close to the source. 

The physical problem can be reduced by pushing the 
outer boundary to larger radius, and the mathematical 
problem will be alleviated at least by making the domain 
of dependence of the initial data larger. There is empiri¬ 
cal evidence that the dominant error in astrophysics-type 
GR simulations is from the boundary, and that it can be 
alleviated for example by using nested boxes mesh refine¬ 
ment to push the outer boundary further out 0,0,0 • 

A radical approach to both the mathematical and the 
physical problems is to represent a complete asymptot¬ 
ically flat spacetime on the numerical domain. The no 
ingoing radiation condition can then be imposed exactly 
at past null infinity and the outgoing radiation can 


be read off exactly at future null infinity . Regularis¬ 
ing the field equations at infinity requires embedding the 
Einstein equations into a much larger system called the 
conformal field equations 0] , but this system has not yet 
been used for astrophysical simulations. Alternatively, 
the Einstein equations can be re gula rised at using 
retarded time as a null coordinate , but this approach 
cannot be used in the central strong field region because 
the null cones form caustics. 

Finally, space can be compactified on standard space¬ 
like slices going out to 7°. If the initial data are sta¬ 
tionary outside a central region one can then get away 
with a pseudo-regularisation of the field equations in 
their discretised form at the compactification boundary 
[1111113 In the continuum equations no gravitational 
waves reach and no evolution takes place there. In 
the discretised equations spurious radiation must be sup¬ 
pressed by artificial dissipation (or the implicit dissipa¬ 
tion of the discretisation) before reaching the compacti¬ 
fication boundary. 

One key ingredient of the conformal field equations 
approach is the use of asymptotically null time slices 
which terminate at , together with a compactifica¬ 
tion of the radial coordinate. It has been suggested to 
use such coordinates in numerical relativity, but to in¬ 
troduce an artificial outer boundary at some very large 
radius and to use the standard form of the Einstein equa¬ 
tions [Ilia. Fig.ui illustrates this approach. Ingo¬ 
ing (backscattered) waves at large radius are increasingly 
lost in numerical simulations in this approach because 
they are not resolved on any finite resolution grid. It 
is sometimes not appreciated that all approaches that 
put at a finite coordinate value, in particular the 
conformal field equations and the Gauchy-characteristic 
matching approaches, suffer from the same shortcom¬ 
ing. For example, in Minkowski spacetime with metric 
ds^ = —dUdV -I- R^dUP' the conformal approach would 
introduce a compactification U = tanu, V = tanu on 
a grid with constant resolution Au and Av. This im¬ 
plies AV = sec^ Au ~ V~‘^Av (as V -i- oo), so that 
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FIG. 1: A schematic spacetime diagram to illustrate that if 
the artificial outer boundary is pushed out to sufficiently large 
radius (here R = 2000M) outgoing gravitational radiation 
can be read off at very large radius (here R = lOOOM) for 
a long time (here t = lOOOM) without being contaminated 
by unphysical BCs at the artificial outer boundary (assuming 
that no information travels faster than light). 


any ingoing waves with period / are no longer resolved 
for V > (/Aw)“^/^. A similar argument holds for the 
characteristic approach on outgoing null cones. 

In this paper, we shall restrict ourselves to the toy 
model of the wave equation on Minkowski or Reissner- 
Nordstrom spacetime. In Sect. HH we investigate in some 
detail the relationship between radial stretching or com- 
pactification, bending up the slices, the coordinate speeds 
of light, and the global structure of the slices. We pro¬ 
pose a one-parameter family of asymptotically null slices 
comprising flat slices at one extreme and hyperboloidal 
slices at the other. 


In Sect. cni we carry out numerical evolutions of an 
outgoing wave packet in the spherical wave equation 
in ADM-like form on the Minkowski and Schwarzschild 
spacetimes. We compare the accuracy and the numeri¬ 
cal cost of evolving on different asymptotically null slices 
with the standard slicing, for realistic values of the nu¬ 
merical parameters. In the Schwarzschild simulations we 
also look for power-law tails and find that they are cor¬ 
rectly resolved. Sect. nvi contains our conclusions. 


II. MATHEMATICAL ANALYSIS 

A. Spherically symmetric slices in Minkowski 
spacetime 

The metric of Minkowski spacetime in standard spher¬ 
ical polar coordinates is 

ds2 = _dr2-bdi?2 + i?2^I]2^ (1) 

where = dO^ + sin^ 9 . 

Consider the flat space wave equation as a toy model 
for the linearised Einstein equations. An outgoing wave 
(j) emanating from a central region can be approximated 
as 

^max ^ 

EE (2) 

l—O m——l 

for some small finite (maxj where the Yim are spherical 
harmonics and U = T — R is retarded time. To resolve 
this, constant spatial resolution Ax, Ay, Az is not re¬ 
quired, but only constant angular resolution A9, A(p. 
Similarly, constant resolution in R and T separately is 
not required, but only constant resolution in U. If the 
slices of constant coordinate time t “bend up” such that 
they approximate slices of constant U asymptotically as 
R^ oo, then the radial grid spacing AR can be allowed 
to diverge in this limit. 

To study this quantitatively, we make the change of 
coordinates 

t = T-FiR), (3) 

R = R{r), (4) 

where F{R) controls the slicing, and R{r) the spatial 
coordinates on the slices. The new slicing is static, that 
is. Lie-dragged by the Killing vector dfdT. 

The metric in 3-1-1 in form in the new coordinates is 

ds^ = —a^dt^ -I- 7rr(dr -|- -I- qegdH^, (5) 

with the lapse, shift, 3-metric and extrinsic curvature 
given by 


a 

1 

(l_ir/2)i’ 

(6) 

Z?” 

F' 

(7) 

'Jrr 

= R''^{1-F'"^), 

(8) 

lee 

= R\ 

(9) 

Krr 

FI'R'2 

(10) 

Kee 

F'R 

(11) 
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Note that a > 1. Here and in the following F' = dF/dR 
and R' = dR/dr. The coordinate light speeds in the 
radial and tangential directions are 


c± 


ce 


dr 
dt 

— - I 
dt V jee 



/?'■ = ± 


R'{1tF')'^ 


( 12 ) 


1 

R' 


(13) 


We can also express the radial coordinate speeds of light 
dr/dt = c± directly in terms of dR/dT = C± in the form 

= R'{CZ^ - F'\ (14) 


and this relation holds for any spherically symmetric met¬ 
ric. 


B. Stability conditions 

We consider a numerical stencil centred at the grid 
point Xq at time tg, and for simplicity of presentation 
we shift the origin of the coordinates so that Xg = 0 
and tg = 0. We model the stencil as the square box 
—Ax* < X* < Ax* for z = 1 , 2 ,3. The time step is At, and 
for simplicity of notation we assume a two-level scheme 
and a three-point stencil in each spatial direction. 

The intersection of the past light cone of the grid point 
at X* = 0 and t = At with the time slice t = 0 is the 
ellipsoid 


The factor 1 / sin 0 in the third condition reflects the 
standard problem with spherical coordinates on the pr¬ 
axis unrelated to our coordinates and which we ignore 
here. 

The CFL condition, which is necessary for stability, 
can be written as 


|c+| < k 


Ar 

'At' 


|c_| < k 


Ar 

At’ 


ct < k 


AO 

'At 


( 20 ) 


for k = 1, where 


_ a 1 

E ^ i?(l-F'2)i/2- 


( 21 ) 


For typical explicit finite differencing methods for lin¬ 
ear hyperbolic equations, a sufficient stability condition 
is then of the same form for some constant k which de¬ 
pends on the equation and its discretisation (allowing for 
any number of time levels and size of stencil). This will 
generalise to any hyperbolic formulation of the Einstein 
equations. Note that ct eg unless a = 1, and that ct 
and not eg is the relevant quantity for stability. 


C. Slicing, stretching, compactification, and 
lightcones 

The assumption that the slice is spacelike is equivalent 
to F' < 1, and the assumption that it becomes asymp¬ 
totically null is equivalent to F' ^ 1. For simplicity let 
us consider functions 1 — F' that decay as a power of R 
at large R, or 


^ijix’' + /3*At)(x'^ -I- (3^ At) = a^At^. (15) 

Note that the origin of this ellipsoid is shifted from x* = 0 
if /?* yf 0 , and its principal axes are not aligned with the 
coordinate axes if 7 ^ is not diagonal. 

The light ellipsoid hts into the stencil box (the 
Courant-Friedrichs-Levy (CFL) condition) if and only if 
in each of the three coordinate directions i the maximal 
value of X* on the ellipsoid is less than Ax* and the min¬ 
imal value is greater than —Ax*. The solution of the 
resulting extremisation problem gives 

— >q;V^+1/3*1 (16) 

for i = 1,2,3 (no summation over i). This result is valid 
also for non-diagonal ■jij. 

Restricted to spherical symmetry, these conditions be¬ 
come 


Ar 



(17) 

At 

> 


A9 


a 

(18) 

'm 

> 

R' 

Aif 

> 

a 

(19) 

At 

i? sin 0 ' 


1 - F'{R) ~ R-** (22) 

as i? ^ 00 , where n > 0. Without loss of generality we 
assume that the numerical grid is equally spaced in r. 
Then the CFL condition in the radial direction requires 
that 

TDTl 

~ ^ (23) 

is bounded above as i? ^ 00 . To preserve accuracy, it 
should also be bounded below. (If c+ ^ 0 as i? ^ 00 , an 
outgoing wave pulse would slow down, become increas¬ 
ingly narrow, and be damped by numerical dissipation.) 
Therefore we require 

C-y ~ ~ const., (24) 

which on integration gives us 

r~i?i-**. (25) 

This means that for 0 < n < 1 we have radial stretching 
i?(r) ~ r ^ 00 , (26) 

and for n > 1 we have radial compactiheation 
i?(r) ~ (/-r)-i/(**-i), r^/_. 


( 27 ) 
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The special case n = 1 gives us 


and the corresponding Lorentz factor is 


Rir) 


oo, 


also a case of stretching. 
We see that 


c_ ~ — 


2R' 


(28) 


(29) 


as i? —*■ oo for either stretching or compactification. This 
means that ingoing wave pulses move very slowly (in r 
and t) at large radius, and that they become wider (in r) 
as they move in. Although they are insufficiently resolved 
at large radius, they do not pose a stability problem for 
the finite differencing scheme. We also see that 

CT ~ (30) 


as i? —> oo, and so with constant angular resolution the 
CFL condition in the tangential directions requires n <2. 

In summary, the ansatz (EH) for the slicing requires 
0 < n < 2 for the coordinate speeds in all directions 
to be bounded. Keeping c+ bounded below as well as 
above means that 0 < n < 1 requires radial stretching 
and 1 < n < 2 requires radial compactification. 

As a matter of convention, we always set R(0) = 0, 
R(0) = 0 and i?'(0) = 1, so that c+(0) = 1. For regularity 
we require R(r) to be odd and F(R) to be even. We 
introduce a length scale L so that 


as R oo. We demand that c+(oo) = 1, which implies 


F' ~ 1 - 



(32) 


1 _ 1 
“ V2 \l) 


(34) 


as R ^ oo. 

We see that our n > 0 slices are asymptotically null 
in the precise sense that their Lorentz factor diverges as 
i? —> oo, but surprisingly this does not imply that they 
reach Consider the null coordinates V = T + R = 
t + F{R) + R and U = T—R = t+F(i?) —i? on Minkowski 
spacetime. Noting that F{R) os R as R ^ oo, we find 
that V os 2R ^ oo as i? ^ oo on a slice of constant t. 
For 17 on a slice of constant t, we find that 


dU 

dR 


= F'{R) -los- 

t=const. 



(35) 


Integrating this relation over R, we find that for n > 1, 
U approaches a finite value as i? —> oo (or equivalently, 
as we have seen, V oo) on a slice of constant t, and so 
the slice reaches a point on For n<l,17^—ooas 
K — > oo on a slice of constant t, and so the slice reaches 
If our slices are plotted on the standard conformal 
diagram, they approach horizontally for n = 0, but 
tangentially to for 0 < n < 1. For 1 < n < 2 they 
intersect 

For comparison we determine the asymptotic be¬ 
haviour of the standard conformal compactification of 
Minkowski M- With U and V as above, and u = 
(t —r)/2 and v = (t-|-r)/2, this is given hy U = tanu and 
V = tanu. (For consistency with the rest of this paper 
u and V have finite range and U and V infinite range). 
Clearly this is is not of the restricted form we have con¬ 
sidered here and which gives a static slicing. But we can 
show that as i? —> oo on a slice of constant t, 


Simple closed form expressions for F{R) and i?(r) with 
the required asymptotic behaviour are given in App. \B 
A typical set of coordinate speeds is plotted in Fig. |21 
In summary, compactifying the radius without bend¬ 
ing up the slices would mean losing the outgoing radi¬ 
ation. Bending up the slices without radial compactifi¬ 
cation would increase the radial outgoing characteristic 
speeds without limit and so violate the CFL condition in 
the radial direction. Bending up the slices without keep¬ 
ing the angular resolution constant at large radius would 
violate the CFL condition in the angular directions. 


dR 


~ 1 - 

t=const. 



(36) 


which we can interpret as n = 2 but where L{t) now 
depends explicitly on the slice. In this compactification 
c± = ±I exactly, and so both ingoing and outgoing waves 
are represented accurately. The lines of constant (r, 9, tp) 
are no longer Killing trajectories, and so we cannot im¬ 
plement this gauge choice with “symmetry-seeking” co¬ 
ordinate conditions |16| , which include almost all popular 
coordinate conditions. 


D. Global structure 


E. Reissner-Nordstrom spacetime 


The exponent n can be given a geometric interpreta¬ 
tion as follows. There are two preferred families of ob¬ 
servers: those normal to the slicing and those along the 
Killing field d/dT = d/dt. The speed of the former rel¬ 
ative to the latter is 



a 


(33) 


To obtain asymptotically null slices in Reissner- 
Nordstrom spacetime, we start from its metric in (for 
example) Kerr-Schild coordinates {T,R), 

ds^ = -{1-f)dT‘^+ 2fdTdR+{l +f)dR^ 

+R‘^ dn^ where /(i?) =‘^ - (^7) 
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Speeds 



FIG. 2: Coordinate speeds c+ and c_ in modified Kerr-Schild 
coordinates, made asymptotically null with n = ^ (solid lines) and 
n = 1 (dashed lines). L = 1 in both cases. 


and make a coordinate transformation of the form (OH) 
to new coordinates {t,r). We choose R{r) to have the 
asymptotic behaviour ED with 0 < n < 2. From CD we 
find that in order to obtain c+(oo) = 1, the asymptotic 
behaviour of F (R) as i? ^ 0 must now be 


. , 4m 8m^ — 2g^ 




(38) 


for 0 < n < 2. For n < 2 the R~^ term in this expression 
can be omitted, and for n < 1 the R~^ term can also 
be omitted, as they give only o(l) contributions to c+. 
Conversely, we see that for n > 1 the slicing must be cor¬ 
rected for the total mass of the isolated central object. 
In a dynamic 3D spacetime this mass would presumably 
be the ADM mass for n = 1 and the Bondi mass for 
n > 1. Although our calculations are limited to spherical 
symmetry, the analogy between charge and angular mo¬ 
mentum suggests that with n = 2 in 3D the slicing would 
have to be corrected for the Bondi angular momentum 
as well. 

Finally we note that (EOl) holds also in Reissner- 
Nordstrom spacetime. More precisely, 


ct ^ 


1 

V2R 



(39) 


for 0 < n < 2 with any m and q, and 


Ct ^ 


y 2 L’ 


L^ = L^ + 24rrF - 


(40) 


for n = 2 as i? —> oo. 


F. Characterisation of n in terms of 3-|-l variables 

We have already characterised the slicing geometrically 
in terms of its Lorentz factor. Another geometric char¬ 
acterisation is in terms of the metric and intrinsic curva¬ 
ture of the slices. Applying the coordinate transforma¬ 
tion (1^ with the asymptotic behaviour mm to the 


metric Ep. we find that for Minkowski with 0 < n < 2 
and Reissner-Nordstrom with 0 < n < 2, 


a 






Kee 

K 


'ygg = R^, 

n /i?\ 

V^\l) 

V2\l) ’ 

n (RY~^ 

2V2L \l) 



(41) 

(42) 

(43) 

(44) 

(45) 


in the limit R —>■ 00 . (For n = 2 with m,q ^ 0, the lead¬ 
ing order coefficients change from the ones given here). 
R is defined geometrically as the area radius. To lead¬ 
ing order in R, there is no effect from the mass m. As 
R 00 , K = K'"r + 2K^g is dominated by K^r for all 
n. For n = 0, K vanishes identically, for 0 < n < 2 it 
goes to zero as R ^ 00 , while for n = 2 it approaches a 
constant. 


III. NUMERICAL EXPERIMENTS 

A. Numerical method 

We solve the wave equation restricted to spherical sym¬ 
metry on Minkowski and Schwarzschild spacetime. Our 
initial data is an outgoing Gaussian pulse near the origin 
(but outside the black hole). On Minkowski spacetime 
the problem has an overall scale-invariance, which we fix 
by setting the width of the pulse to Ai? ~ 1. In the 
Schwarzschild simulations, we also set m = 1. 

In astrophysical simulations we would need to maintain 
a fairly constant resolution in a finite inner region where 
the dynamics take place. As a testbed for 3-dimensional 
simulations, we therefore use a resolution of Ar = 0.1, 
which, with our convention R'(0) = 1, corresponds to 
Ai? ~ 0.1 in the central region, or ^ 10 grid points over 
the width of the pulse. We therefore choose L or i so 
that R'(r) ~ 2 at i? = 10, which means that resolution 
at that radius is still half of that at the centre. 

We solve the wave equation in the continuum form 
(IB3IB4II . We use 2nd and 4th-order accurate finite dif¬ 
ferencing schemes proved to be stable in na , with a 
Courant factor of At/Ar = 0.4. A regular origin is dealt 
with by imposing regularity, and taking appropriate lim¬ 
its in the right hand sides. For black hole excision we use 
outflow BCs. 

We estimate the error and check that the results con¬ 
verge to the expected order by comparing the three 
resolutions Ar = 0.1, 0.05 and 0.025. An example 
is Fig. 121 where we plot the local error e(r, G), where 
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= t{R = Rt.,T = i?*) is the value of t when the maxi¬ 
mum of the Gaussian pulse is at physical radius R — R^, 
in a given slicing. 

B. Outer boundary treatment 

If the outer boundary is at i?max = oo, we impose BCs 
from the exact solution. For the 2nd-order accurate code, 
this simply means imposing (/> = 11 = 0 at the outermost 
grid point i = N where R = oo. In the 4th-order accurate 
code on Minkowski we also impose the exact solution at 
i = N — 1. On Schwarzschild, where no exact solution 
exists, we copy R(j) and i?n from i = N — 2 to i = N— 1. 
(This is only first-order accurate, but in practice does not 
affect results away from the boundary.) 

If the outer boundary is at finite physical radius i?maxj 
we impose maximally dissipative BCs (MDBCs) on (f) and 
n, with lower-order terms representing the l/R falloff. 
This continuum BC perfectly represents outgoing waves 
in spherical symmetry on flat spacetime, but in curved 
spacetime and/or for non-spherical waves gives rise to an 
unavoidable small continuum reflection. The finite dif¬ 
ferencing implementation of this BC also gives rise to 
numerical reflections which converge away with increas¬ 
ing resolution. 

With n > 0 the outer BC is “almost” an outflow 
boundary. A second possible BC is therefore to increase 
the radial shift by a small amount, leaving the other met¬ 
ric coefficients unchanged, to turn the outer boundary 
into a genuine outflow boundary. The background space- 
time is then no longer a vacuum solution of Einstein’s 
equations but the error is small. This is a simple variant 
of Misner’s idea M of introducing a small cosmological 
constant to create a cosmological horizon that represents 
a null outer boundary at constant R. 

As a third BC, we have also tested an outer buffer zone 
in which c+ goes smoothly to zero. This can be achieved 
by reducing n of the slicing and/or increasing n of the 
radial coordinate transformation. In either case the ef¬ 
fect is that outgoing waves become “blue-shifted” with 
respect to the grid spacing and are then suppressed by 
artificial dissipation. We consider the buffer not as part 
of the physical spacetime but rather as a BC. For con¬ 
venience of implementation, with n > 1 in the physical 
region we reduce n for the slicing in the buffer, and with 
n < 1 in the physical region we increase n of the radial 
coordinate transformation (to a value > 1) in the buffer. 
This means that for all values of n in the physical region 
the buffer goes to oo. 

C. Spherical wave equation on Minkowski 

We evolve an outgoing scalar field pulse with different 
slicings with the aim of quantifying the error and estab¬ 
lishing what fraction of grid points can be saved by us¬ 
ing asymptotically null coordinate slices, and what effect 


rescaled error plot for n=1.0 



FIG. 3: The Error in 0 with n = 1 and L = 5.77, at = 10, 
100 and 1000. e(t, r+) at Ar = 0.1 and 16e(r, t+) at Ar = 0.05 
are both plotted, but are indistinguishable here, indicating 4th- 
order convergence. The absolute error e decreases with i?*, but 
the relative error e/R increases. 


they have on the error as the wave goes out to very large 
R. We present results for n = 1/2,1,3/2, 2. The results 
for n < 1 have been obtained with MDBCs at i? > 1000, 
and the results for n > 1 by imposing (/ = 11 = 0 at 
and, in 4th-order accuracy, the exact solution at the grid 
point just inside the boundary. 

We evolve up to the time when the outgoing wave 
reaches the outer boundary. Until then solutions for all 
n converge with the expected order (2nd or 4th), and so 
we have an accurate estimate of the numerical error. The 
number of grid points required to reach a given physical 
radius i?max and the norm of the finite differencing error 
when the wave has gone out to three selected values i?* of 
R are summarised in Table HI Fig. El is a typical example 
of the error plots. 

The larger the value of n, the smaller the number of 
grid points required. The error when the wave reaches 
i? = 10 is similar for all values of n as one would expect 
from the fact that all grids are similar in the central re¬ 
gion R % 10, but by the time the wave reaches R = 1000 
all n > 0 slicings do better than n = 0 , with n = 1 an or¬ 
der of magnitude more accurate than n = 0. We believe 
that this could be explained by the much smaller number 
of time steps required for the wave to reach R — 1000, 
for example 338 for n = 1 compared to 10000 for n = 0. 
This would tend to reduce phase error, and in dissipative 
schemes also amplitude error. 


D. Spherical wave equation on Schwarzschild 

Fields in a curved asymptotically flat spacetime gener- 
ically decay with power law tails because of backscatter 
[T^. In order to see how well our coordinates can cope 
with such small physical effects at large physical radii, we 
evolve the spherical wave equation on the Schwarzschild 
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Slice 

Grid points to cover 
-f^max ~ 1000 -Rmax “ OO 

Time steps to reach 
Rt, = 1000 R* = oo 

|e(',tt)| with 2nd order code at 
R, = 10 77* = 100 77. = 1000 

|e(-,7*)| with 4th order code at 
77. = 10 77. = 100 77. = 1000 

n = 0 
n = 1 

n = 1 

1 

n — 2 

10000 oo 

893 oo 

338 oo 

207 224 

165 167 

25000 oo 

2300 oo 

1000 oo 

624 667 

390 418 

1.8 -lO"® 1.7-10"® 5.0 -lO-"* 
1.5 • 10“® 7.5 • 10"'‘ 2.5 • 10“'‘ 
2.3 • 10“® 8.0 • 10"^ 1.5 • lO-'^ 

2.5-10“® 1.1-10“® 2.5-10“'‘ 
1.5 - 10“® 6.0 - 10“'‘ 5.0 - 10“'‘ 

2.0-10“® 1.9-10“® 1.9-10“® 

1.7- 10“® 1.2 - 10“® 4.5 - 10“® 
3.2-10“® 1.3-10“® 2.0-10“® 
2.3 - 10“® 2.4 - 10“® 4.8 - 10“® 
1.8 - 10“® 1.0 - 10“® 4.0 - 10“® 


TABLE I: Computational costs and errors in evolving a spherical wave on Minkowski spacetime for different slicings. As the 
outgoing wave has amplitude 1/R, a suitable measure of the relative error would be e/R. 


background, on the horizon-penetrating asymptotically 
null slices derived above. A priori it is not clear if tails 
will be represented correctly or not: on the one hand they 
are caused by backscattered ingoing waves, on the other 
hand their frequency in advanced time is low. 

The numerical algorithm is the same as for our 
Minkowski experiments, except for the boundaries. The 
regular origin is replaced by an outflow (excision) bound¬ 
ary inside the black hole. We do not see large errors or 
any sign of instability associated with the excision bound¬ 
ary. We use the 4th-order accurate scheme with the artifi¬ 
cial dissipation parameter tr = 0.007. The results depend 
only weakly on a (including cr = 0). 

At the outer boundary we use three types of BC. The 
first corresponds to what we did in Minkowski spacetime, 
that is, we impose MDBCs at finite i?max or “exact” 
BCs at oo. This gives good tail results if we put the 
outer boundary far enough out. As a second type of 
boundary we have increased the shift near the boundary 
slightly to obtain an outflow BC. As a third type of BC 
we have added a 20% buffer in which c+ goes to zero 
and R reaches infinity (even if n < 1 in the physical 
region). The case n = 0 with such a buffer is similar to 
compactification at used in d, 0,0 • 

As we do not have an exact solution, we have read 
off initial data (j) and 11 on a given initial slice from 
what would be an exact solution on Minkowski where the 
Minkowski coordinates T and R are identified with Kerr- 
Schild coordinates. This means that we have slightly dif¬ 
ferent initial data for different slicings. We have therefore 
used a high resolution simulation on each slicing as a ref¬ 
erence solution to decide when power-law tails are lost. 
We see that at i? = 500 the power is not exactly 3 as 
expected for constant R, but varies between 3 and 2 (the 
power expected at ,%+). 

The errors and computational costs on the domain of 
dependence of the initial data are similar to those given 
in Tabled Therefore we concentrate on power law tails 
at late times. Our results are summarised in Table HD 
To check for tails, we plot 4>{t) for “observers” at i? = 10 
and R = 500. An example is Fig.^ It is worth keeping 
in mind at what time the boundary can have a physi¬ 
cal influence on the observers. For n = 0, an observer 
near the centre leaves the domain of dependence of the 
initial data at t ~ Rmux, while any continuum reflec¬ 
tion of the initial outgoing pulse reaches the observer 



FIG. 4: Log-log plot of <f> versus t at i? = 10 on Schwarzschild, 
using n = 1 slices, with an outer boundary itmax ~ 10®. The 
straight line is with the amplitude but not the power adjusted 
to the data. This is the expected power-law tail. The corresponding 
entry in Table ITTl is in bold face. 


at t ~ 2i?max- With n > 0, one can roughly think of 
the slices as null, and of outgoing waves as going out al¬ 
most instantaneously, so that both events happen close 
together at t ~ 2i?niax- We find that for n = 0, the tail 
is lost at exactly the time when the continuum reflection 
of the outgoing wave reaches the observer (at i? = 10 
and 500 respectively). For n > 0, the tail can be lost 
both earlier (presumably because numerical boundary er¬ 
ror travels faster than light) or later (presumably because 
at very large boundary radius the continuum reflection 
is small and the numerical reflection is also small). 

Our best results are obtained for n = 1 with MDBCs 
imposed at i?niax — 10®? which requires only 722 grid 
points for a central resolution of AR = 0.1. From the 
Minkowski results we see that this choice would also be 
among the best ones for minimising the error in the out¬ 
going wave. 


IV. CONCLUSIONS 

The essential motivation for using either the conformal 
field equations or null coordinates is their ability to sim¬ 
ulate an asymptotically flat spacetime out to on a fi- 
















Slice 

-Rmax 

Grid points 

ti with MDBCs 

or exact BCs at oo 

iZ = 10 7Z = 500 

ti with outflow BCs 

7Z = 10 iZ = 500 

ti with buffer BCs 

iZ = 10 iZ = 500 

n — Q 

1000 

10086 

2020 

1530 


> 2100 

> 2100 

n = 1 

1000 

884 

~ 1700 

- 900 

~ 1800 

~ 900 

> 2100 

> 2100 


10® 

28271 

> 10"* 

> 10^ 

> 2800 

> 2800 



n = 1 

1000 

324 

~ 900 

- 500 

~ 600 

~ 200 

- 900 

~ 500 


10® 

722 

> 10^ 

> 10‘‘ 

> 2400 

> 2400 



n = 1 

1000 

194 

~ 200 

* 

~ 200 

* 

~ 200 

* 


oo 

210 

~ 300 

- 150 





n = 2 

1000 

152 

* 

* 

~ 100 

* 

* 

* 


oo 

153 

~ 300 

- 200 






TABLE II: Time ti when power-law tails are lost for different slicings, different BCs, two outer radii, and two observer locations. 
Note that T = t + const for these observers at fixed R, where the constant depends on R and the slicing. A * means that 
power-law tails were not seen. A blank means that this combination of slice and BC is not sensible. iZmax is the outer radius 
of the physical region (defined arbitrarily as c+ > 1/2), not of the buffer. With buffer BCs the number of grid points increases 
by about 20% from the value stated in the third column. 


nite domain. However, this can be done accurately with 
a finite number of grid points only in situations where 
ingoing radiation, including backscatter, can be progres¬ 
sively neglected at large radius. Conversely, in any phys¬ 
ical situation where these approaches work, one should 
be able to simulate the spacetime almost up to by 
using asymptotically null slices without regularising the 
field equations at infinity min 

From the requirement that the outgoing coordinate 
speed of light should neither go to zero (the code be¬ 
comes inefficient and too dissipative) nor to infinity (ex¬ 
plicit schemes become unstable, and implicit schemes are 
likely to become inaccurate), we have established a quan¬ 
titative relationship between the rate at which the slices 
become null at large radius and the rate at which the 
radial coordinate is compactified. 

This rate can be characterised by a parameter n in 
the range 0 < n < 2, where on the one hand the phys¬ 
ical radius (area radius) R and coordinate radius r are 
asymptotically related by r ^ and on the other 

hand the Lorentz factor of the slices asymptotically di¬ 
verges as r ~ and the extrinsic curvature of the 

slices scales as AT ^ 

The extreme value n = 2 corresponds to the hyper- 
boloidal slices which have traditionally been used both 
in drawing the standard conformal diagram of Minkowski 
or Schwarzschild spacetime, and in numerical evolutions 
using the conformal field equations. For 0 < n < 1 the 
range of the coordinate r is infinite (we prefer to speak of 
stretching rather than of compactification) and surpris¬ 
ingly the slices terminate at For 1 < n < 2 the range 
of r is finite and the slices intersect as expected. 

Our numerical simulations of the spherical wave equa¬ 
tion on Minkowski and Schwarzschild spacetimes have 
confirmed our expectations. For a given (very large) ra¬ 
dius of the outer boundary and a given (realistic) ac¬ 
curacy requirement as the outgoing wave reaches that 


boundary, simulations in our n > 0 coordinates require 
dramatically fewer grid points. To compare two extreme 
examples, evolving an outgoing wave with wavelength 
AR ^ 1 out to R — 1000 with a relative error of 20% 
requires 10000 gridpoints (and 25000 timesteps) in stan¬ 
dard coordinates using 4th-order finite differencing, while 
only 165 grid points (and 390 time steps) with our n = 2 
coordinates give a relative error of only 4%. 

Although n > 0 coordinates do not represent in¬ 
going waves accurately, power-law tails of waves on 
Schwarzschild are correctly represented, until numerical 
error from the boundary becomes bigger than the tail 
signal. Lower values of n represent the tails correctly for 
longer time because numerical error from the boundary 
propagates in less rapidly. This suggests that in appli¬ 
cations to numerical relativity all values of n should be 
considered. From our numerical experiments, the pre¬ 
ferred value of n as a compromise between minimising 
grid size, representing outgoing waves, and representing 
power-law tails is n = 1. It also appears that an n = 2 
slicing would have to be adapted to the angular momen¬ 
tum of the central source, whereas n < 2 slicings are only 
affected by the mass. 

In 3D, our method requires approximately constant 
angular resolution as R oo. A promising technique 
for achieving this while avoiding the axis singularity of 
spherical polar coordinates is the use of multiple coordi¬ 
nate patches [UllSilllil. 

In summary, our proposed approach for pushing the 
BCs to very large radius requires only a modification of 
the initial data and the gauge conditions but uses stan¬ 
dard 3-1-1 field equations. It is simpler than the conformal 
field equation method and Cauchy-characteristic match¬ 
ing methods. Although it cannot reach , it can rep¬ 
resent outgoing waves out to arbitrarily large radius on a 
finite grid. Our results also suggest certain improvements 
to existing numerical relativity codes: 
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a) A lower value of n, perhaps n = 1, may be better 
than n = 2 for representing physically interesting scenar¬ 
ios accurately. The conformal and null methods require 
n = 2 because they require a regular conformal spacetime 
metric, and this requires the conformal factor ^ R~^ as 

is approached. Keeping this asymptotic behaviour, 
one may be able to improve their accuracy at large R by 
using a smaller value of n out to some very large value 
of R and then switching to n = 2 to reach . 

b) In 3D implementations of compactification at 
[M113 7 Cartesian spatial coordinates are compactified as 
A* = tana:®, which corresponds to n = 2 compactifica¬ 
tion in our notation. Our results suggest a modification 
of this approach in which the slices are made asymptot¬ 
ically null up to large finite radius, perhaps best with 
n = 1. The resulting code would represent outgoing ra¬ 
diation correctly out to this large finite radius. Outside 
this region the slice could still reach for example by 
switching to n = 2 (the traditional value). 

c) Fixed mesh refinement [1,0 0] ? in 3D typically using 
Cartesian grids, is also widely used to push the outer 
boundary further out. If this is done with nested boxes 
with equal (or comparable) numbers of grid points in 
each box, one automatically has approximately constant 
angular resolution, and radial resolution Aii R, which 
corresponds to n = 1 stretching in our notation. Our 
results suggest combining the nested boxes approach with 
n = 1 slicing (and a global time step), so that outgoing 
waves are resolved all the way out to the outer boundary, 
which can then be put at arbitrarily large radius. 


for 0 < n < 2 with n and 

Gi(i?, L) = — In -1-—^ (A5) 

for the special value n = 1. They obey 

G'^{R,L)=(^ +0{R—^). (A6) 

for all n including n = 1. In particular, 

GUi?,L) = l-^ + 0(i?-4), (A7) 

and the R~^ term in Gg has to be taken into account in 
the n = 2 slicing. It is straightforward to construct F{R) 
from these blocks. We only give an example here. For 
0 < n < 2 with m = q = 0, we can use 

L 

G„(i?,L„). (A8) 

Note that Lq and L„ can differ from each other and from 
L (the constant introduced in R{r)) as long as the corre¬ 
sponding amplitude is adjusted. For n > 1 with m > 0 
and again for n = 2, additional terms have to be added. 
Similarly, to achieve a switchover from n to noo < n at 
a given large R we only need to add G„^ with suitable 
amplitude and scale. 


F{R) = Go{R, Lo)- 


APPENDIX A: CLOSED FORM COORDINATE 
TRANSFORMATIONS 

We give closed-form functions R{r) and F{R) which 
have the required asymptotic behaviour ED and (ISSli . 
For the radial stretch we choose 

/ 1 \i/" 

i!M=r(l + ^) . (=(—) L(A1) 

for the generic case 0 < n < 1, and 

i?(r) =Lsinh— (A 2) 

Lj 

for n = 1. The range of r is 0 < r < oo. For the 
compactification we choose 

for 1 < n < 2. The range of r is now 0 < r < 1. 

The expression for F(R) can be assembled from the 
functions 


APPENDIX B: WAVE EQUATION AND 
BOUNDARY CONDITIONS 

We write the massless wave equation 

VaV^c/> = 0 (Bl) 

in a 3 -I-1 form similar to the ADM form of the Einstein 
equations by defining 

n = -£r,c/> = -n^Vn(/>, (B2) 

where Ua is the future-pointing unit normal on the slices 
of constant t. With F = {d/dt)°' = an°' + /3“ this gives 

Lt4> = - an, (B3) 

An = Cf}n-aDaD'^<f)-aa!^Db(f) + aKn. (B4) 

Here Da is the covariant derivative associated with the 
three-metric ^ab on each slice, and = D^\n.a is the 
acceleration of the observers. 

At the outer boundary the MDBC corresponding to no 
incoming radiation |23| would be, 

n - m'^Daft) = 0, (B5) 


rn 1-n 

Gn{R,L) = - - (i?2 + A) 

1 — n ^ ^ 


L 


(A4) 


where is the outward-pointing unit vector normal to 
the boundary within each constant t slice. In Minkowski 


1 — n 
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this is equivalent to 0 = f{R—T). By contrast an outgo¬ 
ing spherical wave in Minkowski is </> = f{R — T)/R. We 
can achieve this boundary condition by adding a lower- 
order term to the MDBC, namely 


n - m°^Da(j) + Q(j) = Q, 


Q = 


« - /3’~a/w 

\/W 


(B6) 


Note that a MDBC is defined to lead to a non-increasing 
energy in the frozen coefficient, principal part only ap¬ 
proximation, so that this modified BC is still a MDBC. 

Following Ea , we discretise this to 2nd-order accuracy 
as 


= 0, (Bll) 

h'^DXliN+2 = 0, (B12) 

where is defined by 

(B13) 

At the black hole-excision boundary, or where the outer 
boundary is null or spacelike, no BCs are required for 
the continuum equations. As numerical BCs we then use 
extrapolation for all ghost points including 4>n+i^ that is 


IItv — '/¥^Do4’n + QN(t>N = 0, (B7) 

h^DtUN+i = 0, (B8) 

where N is the grid point at r = r^ax and iV -|- 1 is a 
ghost point. For 4-th order accuracy we use 

+ QN(j)N = 0, (B9) 

h^Dl(l}N+2 = 0, (BIO) 


h^Dlcl)N+i = 0 (B14) 

in the 4-th order accurate case, and in the 2nd-order 
accurate case. 

Our artificial dissipation operators are —ah^{D^D-)^ 
at 2nd-order accuracy and ah^{D+D-)^ at 4th order. 
We use extrapolation to populate the additional ghost 
points required at boundaries. 
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